20225421 권다희 HW4

0. Load Image with dcraw

To load linear raw image, we need to convert them into readable linear 16-bit image.
dc = readraw;
readraw: installed ReadRaw as loader for RAW camera images.
im1_r = imread(dc, '../exposure_stack/exposure1.nef', '-w -T -6 -q 3 -o 1');
im2_r = imread(dc, '../exposure_stack/exposure2.nef', '-w -T -6 -q 3 -o 1');
im3_r = imread(dc, '../exposure_stack/exposure3.nef', '-w -T -6 -q 3 -o 1');
im4_r = imread(dc, '../exposure_stack/exposure4.nef', '-w -T -6 -q 3 -o 1');
im5_r = imread(dc, '../exposure_stack/exposure5.nef', '-w -T -6 -q 3 -o 1');
im6_r = imread(dc, '../exposure_stack/exposure6.nef', '-w -T -6 -q 3 -o 1');
im7_r = imread(dc, '../exposure_stack/exposure7.nef', '-w -T -6 -q 3 -o 1');
im8_r = imread(dc, '../exposure_stack/exposure8.nef', '-w -T -6 -q 3 -o 1');
im9_r = imread(dc, '../exposure_stack/exposure9.nef', '-w -T -6 -q 3 -o 1');
im10_r = imread(dc, '../exposure_stack/exposure10.nef', '-w -T -6 -q 3 -o 1');
im11_r = imread(dc, '../exposure_stack/exposure11.nef', '-w -T -6 -q 3 -o 1');
im12_r = imread(dc, '../exposure_stack/exposure12.nef', '-w -T -6 -q 3 -o 1');
im13_r = imread(dc, '../exposure_stack/exposure13.nef', '-w -T -6 -q 3 -o 1');
im14_r = imread(dc, '../exposure_stack/exposure14.nef', '-w -T -6 -q 3 -o 1');
im15_r = imread(dc, '../exposure_stack/exposure15.nef', '-w -T -6 -q 3 -o 1');
im16_r = imread(dc, '../exposure_stack/exposure16.nef', '-w -T -6 -q 3 -o 1');
imgs_r = {im1_r,im2_r,im3_r,im4_r,im5_r,im6_r,im7_r, im8_r,im9_r,im10_r,im11_r,im12_r,im13_r,im14_r,im15_r,im16_r};

0. Load non-linear rendered images

I also load non-linear images (.jpg files)
im1 = imread('../exposure_stack/exposure1.jpg');
im2 = imread('../exposure_stack/exposure2.jpg');
im3 = imread('../exposure_stack/exposure3.jpg');
im4 = imread('../exposure_stack/exposure4.jpg');
im5 = imread('../exposure_stack/exposure5.jpg');
im6 = imread('../exposure_stack/exposure6.jpg');
im7 = imread('../exposure_stack/exposure7.jpg');
im8 = imread('../exposure_stack/exposure8.jpg');
im9 = imread('../exposure_stack/exposure9.jpg');
im10 = imread('../exposure_stack/exposure10.jpg');
im11 = imread('../exposure_stack/exposure11.jpg');
im12 = imread('../exposure_stack/exposure12.jpg');
im13 = imread('../exposure_stack/exposure13.jpg');
im14 = imread('../exposure_stack/exposure14.jpg');
im15 = imread('../exposure_stack/exposure15.jpg');
im16 = imread('../exposure_stack/exposure16.jpg');
imgs = {im1,im2,im3,im4,im5,im6,im7,im8,im9,im10,im11,im12,im13,im14,im15,im16};

1.1 Linearize - with uniform weight

For High Dynamic Range imaging, we need linear images. Therefore, I linearize the .jpg images. For linearization, I solve the linear equation introduced in the lecture.
w = weighten('unif','jpg');
B = log([1/2048,1/1024,1/512,1/256,1/128,1/64,1/32,1/16,1/8,1/4,1/2,1,2,4,8,16]);
l = 150;
g = zeros(256,3);
lE = zeros(300,3);
 
for c = 1:3
Z = image_stack(imgs,c,300);
[g_c, lE_c] = gsolve(Z,B,l,w);
g(:,c) = g_c;
lE(:,c) = lE_c;
end
plot(g(:,1),1:256);
title('The curvature of function g with uniform weight');
xlabel('ln exposure');
ylabel('pixel value');

- Convert non-linear images into linear ones

% for linear domain
lin_imgs = cell(1,16);
 
 
for c = 1:3
for k = 1:size(imgs,2)
g_tmp = g(:,c);
lin_imgs{k}(:,:,c) = exp(g_tmp(imgs{k}(:,:,c)+1));
end
c
end
c = 1
c = 2
c = 3

1.1 Linearize - with tent weight

w = weighten('tent','jpg');
B = log([1/2048,1/1024,1/512,1/256,1/128,1/64,1/32,1/16,1/8,1/4,1/2,1,2,4,8,16]);
l = 150;
g = zeros(256,3);
lE = zeros(300,3);
 
for c = 1:3
Z = image_stack(imgs,c,300);
[g_c, lE_c] = gsolve(Z,B,l,w);
g(:,c) = g_c;
lE(:,c) = lE_c;
end
plot(g(:,1),1:256);
title('The curvature of function g with tent weight');
xlabel('ln exposure');
ylabel('pixel value');

- Convert non-linear images into linear ones

% for linear domain
lin_imgs_t = cell(1,16);
 
for c = 1:3
for k = 1:size(imgs,2)
g_tmp = g(:,c);
lin_imgs_t{k}(:,:,c) = exp(g_tmp(imgs{k}(:,:,c)+1));
end
c
end
c = 1
c = 2
c = 3

1.1 Linearize - with Gaussian filter

x = (0:255)/255;
w = exp(-4.*(x-0.5).^2/(0.5.^2));
B = log([1/2048,1/1024,1/512,1/256,1/128,1/64,1/32,1/16,1/8,1/4,1/2,1,2,4,8,16]);
l = 150;
g = zeros(256,3);
lE = zeros(300,3);
 
for c = 1:3
Z = image_stack(imgs,c,300);
[g_c, lE_c] = gsolve(Z,B,l,w);
g(:,c) = g_c;
lE(:,c) = lE_c;
end
plot(g(:,1),1:256);
title('The curvature of function g with gaussain weight');
xlabel('ln exposure');
ylabel('pixel value');

- Convert non-linear images into linear ones

% for linear domain
lin_imgs_g = cell(1,16);
 
for c = 1:3
for k = 1:size(imgs,2)
g_tmp = g(:,c);
lin_imgs_g{k}(:,:,c) = exp(g_tmp(imgs{k}(:,:,c)+1));
end
c
end
c = 1
c = 2
c = 3

Merge Exposure Stack into HDR Image

With the linearized images, I made HDR images. There are 3 options for merging exposure stack.

1. uniform weight filter, rendered images with linear domain

w = weighten('unif','jpg');
unif_lin_rendered = zeros(size(im1,1),size(im1,2),size(im1,3));
 
for c = 1:3
nom = zeros(size(im1,1),size(im1,2));
denom = zeros(size(im1,1),size(im1,2));
for k = 1:size(imgs,2)
nom = nom + w(imgs{k}(:,:,c)+1).*im2double(lin_imgs{k}(:,:,c))./exp(B(k));
denom = denom + w(imgs{k}(:,:,c)+1);
end
unif_lin_rendered(:,:,c) = im2double(nom./denom);
end
hdrwrite(unif_lin_rendered,'../results/unif_lin_rendered.hdr');

2. uniform weight filter, rendered images with log domain

w = weighten('unif','jpg');
unif_log_rendered = zeros(size(im1,1),size(im1,2),size(im1,3));
 
for c = 1:3
nom = zeros(size(imgs{k},1),size(imgs{k},2));
denom = zeros(size(imgs{k},1),size(imgs{k},2));
for k = 1:size(imgs,2)
nom = nom + w(imgs{k}(:,:,c)+1).*(log(im2double(lin_imgs{k}(:,:,c))) - B(k));
denom = denom + w(imgs{k}(:,:,c)+1);
end
unif_log_rendered(:,:,c) = exp(nom./denom);
end
hdrwrite(unif_log_rendered,'../results/unif_log_rendered.hdr');

3. uniform weight filter, raw images with linear domain

w = weighten('unif','raw');
unif_lin_raw = zeros(size(im1_r,1),size(im1_r,2),size(im1_r,3));
 
for c = 1:3
nom = zeros(size(imgs_r{k},1),size(imgs_r{k},2));
denom = zeros(size(imgs_r{k},1),size(imgs_r{k},2));
for k = 1:size(imgs_r,2)
nom = nom + w(imgs_r{k}(:,:,c)+1).*im2double(imgs_r{k}(:,:,c))./exp(B(k));
denom = denom + w(imgs_r{k}(:,:,c)+1);
end
unif_lin_raw(:,:,c) = im2double(nom./denom);
end
hdrwrite(unif_lin_raw,'../results/unif_lin_raw.hdr');

4. uniform weight filter, raw images with log domain

w = weighten('unif','raw');
unif_log_raw = zeros(size(im1_r,1),size(im1_r,2),size(im1_r,3));
 
for c = 1:3
nom = zeros(size(imgs_r{k},1),size(imgs_r{k},2));
denom = zeros(size(imgs_r{k},1),size(imgs_r{k},2));
for k = 1:size(imgs_r,2)
nom = nom + w(imgs_r{k}(:,:,c)+1).*(log(im2double(imgs_r{k}(:,:,c))) - B(k));
denom = denom + w(imgs_r{k}(:,:,c)+1);
end
unif_log_raw(:,:,c) = exp(nom./denom);
end
hdrwrite(unif_log_raw,'../results/unif_log_raw.hdr');

5. tent weight filter, rendered images with linear domain

w = weighten('tent','jpg');
tent_lin_rendered = zeros(size(im1,1),size(im1,2),size(im1,3));
 
for c = 1:3
nom = zeros(size(im1,1),size(im1,2));
denom = zeros(size(im1,1),size(im1,2));
for k = 1:size(imgs,2)
nom = nom + w(imgs{k}(:,:,c)+1).*im2double(lin_imgs_t{k}(:,:,c))./exp(B(k));
denom = denom + w(imgs{k}(:,:,c)+1);
end
tent_lin_rendered(:,:,c) = im2double(nom./denom);
end
hdrwrite(tent_lin_rendered,'../results/tent_lin_rendered.hdr');

6. tent weight filter, rendered images with log domain

w = weighten('tent','jpg');
tent_log_rendered = zeros(size(im1,1),size(im1,2),size(im1,3));
B = log([1/2048,1/1024,1/512,1/256,1/128,1/64,1/32,1/16,1/8,1/4,1/2,1,2,4,8,16]);
 
for c = 1:3
nom = zeros(size(imgs{k},1),size(imgs{k},2));
denom = zeros(size(imgs{k},1),size(imgs{k},2));
for k = 1:size(imgs,2)
nom = nom + w(imgs{k}(:,:,c)+1).*(log(im2double(lin_imgs_t{k}(:,:,c))) - B(k));
denom = denom + w(imgs{k}(:,:,c)+1);
end
tent_log_rendered(:,:,c) = exp(nom./denom);
end
hdrwrite(tent_log_rendered,'../results/tent_log_rendered.hdr');

7. tent weight filter, raw images with linear domain

w = weighten('tent','raw');
tent_lin_raw = zeros(size(im1_r,1),size(im1_r,2),size(im1_r,3));
 
for c = 1:3
nom = zeros(size(imgs_r{k},1),size(imgs_r{k},2));
denom = zeros(size(imgs_r{k},1),size(imgs_r{k},2));
for k = 1:size(imgs_r,2)
nom = nom + w(imgs_r{k}(:,:,c)+1).*im2double(imgs_r{k}(:,:,c))./exp(B(k));
denom = denom + w(imgs_r{k}(:,:,c)+1);
end
tent_lin_raw(:,:,c) = im2double(nom./denom);
end
hdrwrite(tent_lin_raw,'../results/tent_lin_raw.hdr');

8. tent weight filter, raw images with log domain

w = weighten('tent','raw');
tent_log_raw = zeros(size(im1_r,1),size(im1_r,2),size(im1_r,3));
 
for c = 1:3
nom = zeros(size(imgs_r{k},1),size(imgs_r{k},2));
denom = zeros(size(imgs_r{k},1),size(imgs_r{k},2));
for k = 1:size(imgs_r,2)
nom = nom + w(imgs_r{k}(:,:,c)+1).*(log(im2double(imgs_r{k}(:,:,c))) - B(k));
denom = denom + w(imgs_r{k}(:,:,c)+1);
end
tent_log_raw(:,:,c) = exp(nom./denom);
end
hdrwrite(im2double(tent_log_raw),'../results/tent_log_raw.hdr');

9. gaussian weight filter, rendered images with linear domain

x = (0:255)/255;
w = exp(-4.*(x-0.5).^2/(0.5.^2));
 
gauss_lin_rendered = zeros(size(im1,1),size(im1,2),size(im1,3));
 
for c = 1:3
nom = zeros(size(im1,1),size(im1,2));
denom = zeros(size(im1,1),size(im1,2));
for k = 1:size(imgs,2)
nom = nom + w(imgs{k}(:,:,c)+1).*im2double(lin_imgs_g{k}(:,:,c))./exp(B(k));
denom = denom + w(imgs{k}(:,:,c)+1);
end
gauss_lin_rendered(:,:,c) = im2double(nom./denom);
end
hdrwrite(gauss_lin_rendered,'../results/gauss_lin_rendered.hdr');

10. gaussian weight filter, rendered images with log domain

x = (0:255)/255;
w = exp(-4.*(x-0.5).^2/(0.5.^2));
 
gauss_log_rendered = zeros(size(im1,1),size(im1,2),size(im1,3));
 
for c = 1:3
nom = zeros(size(imgs{k},1),size(imgs{k},2));
denom = zeros(size(imgs{k},1),size(imgs{k},2));
for k = 1:size(imgs,2)
nom = nom + w(imgs{k}(:,:,c)+1).*(log(im2double(lin_imgs_g{k}(:,:,c))) - B(k));
denom = denom + w(imgs{k}(:,:,c)+1);
end
gauss_log_rendered(:,:,c) = exp(nom./denom);
end
 
hdrwrite(gauss_log_rendered,'../results/gauss_log_rendered.hdr');

11. gaussian weight filter, raw images with linear domain

x = (0:65535)/65535;
w = exp(-4.*(x-0.5).^2/(0.5.^2));
 
gauss_lin_raw = zeros(size(im1_r,1),size(im1_r,2),size(im1_r,3));
 
for c = 1:3
nom = zeros(size(imgs_r{k},1),size(imgs_r{k},2));
denom = zeros(size(imgs_r{k},1),size(imgs_r{k},2));
for k = 1:size(imgs_r,2)
nom = nom + w(imgs_r{k}(:,:,c)+1).*im2double(imgs_r{k}(:,:,c))./exp(B(k));
denom = denom + w(imgs_r{k}(:,:,c)+1);
end
gauss_lin_raw(:,:,c) = im2double(nom./denom);
end
 
hdrwrite(gauss_lin_raw,'../results/gauss_lin_raw.hdr');

12. gaussian weight filter, raw images with log domain

x = (0:65535)/65535;
w = exp(-4.*(x-0.5).^2/(0.5.^2));
 
gauss_log_raw = zeros(size(im1_r,1),size(im1_r,2),size(im1_r,3));
 
for c = 1:3
nom = zeros(size(imgs_r{k},1),size(imgs_r{k},2));
denom = zeros(size(imgs_r{k},1),size(imgs_r{k},2));
for k = 1:size(imgs_r,2)
nom = nom + w(imgs_r{k}(:,:,c)+1).*(log(im2double(imgs_r{k}(:,:,c))) - B(k));
denom = denom + w(imgs_r{k}(:,:,c)+1);
end
gauss_log_raw(:,:,c) = exp(nom./denom);
end
 
hdrwrite(gauss_log_raw,'../results/gauss_log_raw.hdr');

1.3 Evaluation

imshow(im16);
[x1,y1]= ginput(4);
[x2,y2]= ginput(4);
[x3,y3]= ginput(4);
[x4,y4]= ginput(4);
[x5,y5]= ginput(4);
[x6,y6]= ginput(4);
x = [x1,x2,x3,x4,x5,x6];
y = [y1,y2,y3,y4,y5,y6];
patch_avg = zeros(6,1);
patch_actual_avg = zeros(6,1);
 
for i = 1:6
min_x = min(round(x(:,i))); max_x = max(round(x(:,i))); min_y = min(round(y(:,i))); max_y = max(round(y(:,i)));
 
patch = tent_lin_raw(min_y:max_y, min_x:max_x,:);
patch_actual = im16_r(min_y:max_y, min_x:max_x,:);
 
patch_xyz = rgb2xyz(patch,'ColorSpace','linear-rgb');
patch_xyz_actual = rgb2xyz(patch_actual,'ColorSpace','linear-rgb');
patch_avg(i) = mean(patch_xyz(:,:,2),'all');
patch_actual_avg(i) = mean(patch_xyz_actual(:,:,2),'all');
end
lm_model = fitlm(log(patch_avg),log(patch_actual_avg));
plot(lm_model)
title('Linear regression to the log(avg luminance)')
SSE = sum(lm_model.Residuals{:,1}.^2)
SSE = 0.0115
patch_avg_tliraw = patch_avg;
figure, hold on
p1 = plot(log(patch_avg_tlor));
p2 = plot(log(patch_avg_tlir));
p3 = plot(log(patch_avg_ulir));
p4 = plot(log(patch_avg_ulor));
p5 = plot(log(patch_avg_tloraw));
p6 = plot(log(patch_avg_tliraw));
p7 = plot(log(patch_avg_uliraw));
p8 = plot(log(patch_avg_uloraw));
lgd = legend([p1,p4,p6,p8],{'tent/log/rendered','uniform/log/rendered','tent/linear/raw','uniform/log/raw'});
title('Logarithmic plot of average luminances for various HDR imgs')
hold off

By considering these overall evaluations, I will use the "tent weight/logarithm merging/rendered image" for the following tonemapping.

2.1 Photographic Tonemapping

I do tonemapping on both the RGB color space and the XYZ luminance space (only Y).

- For color channel

B = 0.5;
K = 0.8;
I_hdr = zeros(size(im1,1),size(im1,2),size(im1,3));
I_tonemap = zeros(size(im1,1),size(im1,2),size(im1,3));
 
for c = 1:3
I_m = exp(mean(log(tent_log_rendered(:,:,c)+0.00001),'all'));
I_hdr(:,:,c) = K/I_m*tent_log_rendered(:,:,c);
I_white = B * max(I_hdr(:,:,c),[],'all');
 
I_tonemap(:,:,c) = I_hdr(:,:,c).*(1+I_hdr(:,:,c)/(I_white.^2))./(1+I_hdr(:,:,c));
end
I_tonemap = im2uint8(I_tonemap);
imwrite(I_tonemap,'../results/I_tonemap_rgb.jpg');

- For luminance channel Y

xyz = rgb2xyz(tent_log_rendered,'Colorspace','linear-rgb');
xyY = zeros(size(im1,1),size(im1,2),size(im1,3));
xyY(:,:,1) = xyz(:,:,1)./(xyz(:,:,1)+xyz(:,:,2)+xyz(:,:,3));
xyY(:,:,2) = xyz(:,:,2)./(xyz(:,:,1)+xyz(:,:,2)+xyz(:,:,3));
xyY(:,:,3) = xyz(:,:,2);
B = 0.5;
K = 0.8;
I_hdr = zeros(size(im1,1),size(im1,2));
 
 
I_m = exp(mean(log(xyY(:,:,3)+0.00001),'all'));
I_hdr = K/I_m*xyY(:,:,3);
I_white = B * max(I_hdr,[],'all');
 
I_tonemap = I_hdr.*(1+I_hdr/(I_white.^2))./(1+I_hdr);
 
I_tonemap_xyz = zeros(size(im1,1),size(im1,2),size(im1,3));
I_tonemap_xyz(:,:,1) = xyY(:,:,1).*I_tonemap./xyY(:,:,2);
I_tonemap_xyz(:,:,2) = I_tonemap;
I_tonemap_xyz(:,:,3) = I_tonemap.*(1-xyY(:,:,1)-xyY(:,:,2))./xyY(:,:,2);
 
I_tonemap_rgb = xyz2rgb(I_tonemap_xyz);
I_tonemap = im2uint8(I_tonemap_rgb);
imwrite(I_tonemap,'../results/I_tonemap_l.jpg');
The results of the color space are generally blue, and the results of luminescence space are yellow. The larger the K value, the brighter it becomes (the better for this image), and the larger the B value, the greater the contrast (the smaller the influence of B in this image).

With the Photographic Tonemapping algorithm, I like "K=0.8/B=0.5 on RGB space" the most.

2.2 Tonemapping using bilateral filtering

- For color channel

Lij = log(tent_log_rendered);
 
%bilateral filtering
B_filter = zeros(size(im2,1),size(im2,2),size(im2,3));
Lij_nom = zeros(size(im2,1),size(im2,2),size(im2,3));
for c = 1:3
Lij_nom(:,:,c) = rescale(Lij(:,:,c));
end
sigma_r = 0.4;
w = 10;
B_filter = bfilter2(Lij_nom,w,sigma_r);
B_filter_denom = B_filter.*(max(B_filter,[],'all')-min(B_filter,[],'all'))+min(B_filter,[],'all');
Dij = Lij - B_filter_denom;
B_filter2 = 2 .* (B_filter_denom - max(B_filter_denom,[],'all'));
I_tonemap_filter = im2uint8(exp(B_filter2 + Dij));
 
imwrite(I_tonemap_filter,'../results/I_tonemap_filter_c_2.jpg');

- For the luminance

xyz = rgb2xyz(tent_log_rendered,'Colorspace','linear-rgb');
xyY = zeros(size(im1,1),size(im1,2),size(im1,3));
xyY(:,:,1) = xyz(:,:,1)./(xyz(:,:,1)+xyz(:,:,2)+xyz(:,:,3));
xyY(:,:,2) = xyz(:,:,2)./(xyz(:,:,1)+xyz(:,:,2)+xyz(:,:,3));
xyY(:,:,3) = xyz(:,:,2);
Lij = log(xyY);
 
%bilateral filtering
B_filter = zeros(size(im2,1),size(im2,2));
Lij_nom = zeros(size(im2,1),size(im2,2),size(im2,3));
for c = 1:3
Lij_nom(:,:,c) = rescale(Lij(:,:,c));
end
sigma_r = 0.4;
w = 20;
B_filter = bfilter2(Lij_nom(:,:,3),w,sigma_r);
B_filter_denom = B_filter.*(max(B_filter,[],'all')-min(B_filter,[],'all'))+min(B_filter,[],'all');
Dij = Lij(:,:,3) - B_filter_denom;
B_filter2 = 2 .* (B_filter_denom - max(B_filter_denom,[],'all'));
I_tonemap_filter = exp(B_filter2 + Dij);
 
I_tonemap_xyz = zeros(size(im1,1),size(im1,2),size(im1,3));
I_tonemap_xyz(:,:,1) = xyY(:,:,1).*I_tonemap_filter./xyY(:,:,2);
I_tonemap_xyz(:,:,2) = I_tonemap_filter;
I_tonemap_xyz(:,:,3) = I_tonemap_filter.*(1-xyY(:,:,1)-xyY(:,:,2))./xyY(:,:,2);
 
I_tonemap_rgb = im2uint8(xyz2rgb(I_tonemap_xyz));
 
imwrite(I_tonemap_rgb,'../results/I_tonemap_filter_l_2.jpg');
The color space has a strong contrast overall, so the details of the dark part are lost, and the luminance space is not. The image contrast varies depending on the S value, and the result is the best when S=2. However, considering that the result saved as an hdr image and the result converted to uint8 image are quite different, scaling seems to have not been done well.

In conclusion, luminance space with S=2 has the best result.

Complare with the Photographic Tonemapping algorithm, I like the tonemapping with bilateral filter on luminance space (S=2) the most. It is because they quite well preserve the color scale, while maintains the details.

3. Bonus: Implement a different gradient-domain processing

1. Convert to the XYZ space and extract the Y luminance values.

xyz = rgb2xyz(tent_log_rendered,'Colorspace','linear-rgb');
H = imresize(log(xyz(:,:,2)),0.1);

2. Construct the Gaussian pyramid

H_pyramid = genpyramid(H,3);

3. Compute the gradients on each levels of pyramid

H_pyramid_grad_x = cell(1,4);
H_pyramid_grad_y = cell(1,4);
for k = 1:4
H_pyramid_grad_x{k} = zeros(size(H_pyramid{k},1)-1,size(H_pyramid{k},2));
for i = 1:size(H_pyramid{k},1)-1
for j = 1:size(H_pyramid{k},2)
H_pyramid_grad_x{k}(i,j) = H_pyramid{k}(i+1,j)-H_pyramid{k}(i,j);
end
end
end
 
for k = 1:4
H_pyramid_grad_y{k} = zeros(size(H_pyramid{k},1),size(H_pyramid{k},2)-1);
for i = 1:size(H_pyramid{k},1)
for j = 1:size(H_pyramid{k},2)-1
H_pyramid_grad_y{k}(i,j) = H_pyramid{k}(i,j+1)-H_pyramid{k}(i,j);
end
end
end

4. Compute the attenuation function Phi iteratively

beta = 0.8;
Pi_x = cell(1,4);
Pi_y = cell(1,4);
 
%iteratively find Pi
for i = 1:4
alpha_x = mean(sqrt((H_pyramid_grad_x{5-i}).^2),'all').*0.1;
alpha_y = mean(sqrt((H_pyramid_grad_y{5-i}).^2),'all').*0.1;
factor_x = (alpha_x ./sqrt((H_pyramid_grad_x{5-i}).^2)) .* (sqrt((H_pyramid_grad_x{5-i}).^2)/alpha_x).^beta;
factor_y = (alpha_y ./sqrt((H_pyramid_grad_y{5-i}).^2)) .* (sqrt((H_pyramid_grad_y{5-i}).^2)/alpha_y).^beta;
if i == 1
Pi_x{5-i} = factor_x;
Pi_y{5-i} = factor_y;
else
Pi_x{5-i} = imresize(Pi_x{6-i},[size(factor_x,1),size(factor_x,2)],'bilinear').* factor_x;
Pi_y{5-i} = imresize(Pi_y{6-i},[size(factor_y,1),size(factor_y,2)],'bilinear').*factor_y;
end
end

5. Compute the G which represents the attenuated derivates

G_x = H_pyramid_grad_x{1}.*Pi_x{1};
G_y = H_pyramid_grad_y{1}.*Pi_y{1};

6. Solving Poisson equation to get reconstruced pixel intensity.

[imh, imw] = size(H);
im2var = zeros(imh,imw);
im2var(1:imh*imw) = 1:imh*imw;
 
nrow = imh*imw;
ncol = imh*imw;
A = sparse([], [], [], ncol, nrow, 0);
b = sparse(ncol,1);
e = 0;
 
for i = 2:imh-1
for j = 2:imw-1
e = e+1;
A(e,im2var(i,j)) = -4;
A(e,im2var(i-1,j)) = 1;
A(e,im2var(i,j-1)) = 1;
A(e,im2var(i+1,j)) = 1;
A(e,im2var(i,j+1)) = 1;
b(e) = G_x(i,j) + G_y(i,j) - G_x(i-1,j) - G_y(i,j-1);
end
end
for j = 2:imw-1
i = 1;
e = e+1;
A(e,im2var(i,j)) = -4;
A(e,im2var(i,j-1)) = 1;
A(e,im2var(i+1,j)) = 1;
A(e,im2var(i,j+1)) = 1;
 
b(e) = G_x(i,j) + G_y(i,j) - G_y(i,j-1);
end
 
for j = 2:imw-1
i = imh;
e = e+1;
A(e,im2var(i,j)) = -4;
A(e,im2var(i-1,j)) = 1;
A(e,im2var(i,j-1)) = 1;
A(e,im2var(i,j+1)) = 1;
 
b(e) = G_y(i,j) - G_x(i-1,j) - G_y(i,j-1);
end
 
for i = 2:imh-1
j = 1;
e = e+1;
A(e,im2var(i,j)) = -4;
A(e,im2var(i-1,j)) = 1;
A(e,im2var(i+1,j)) = 1;
A(e,im2var(i,j+1)) = 1;
 
b(e) = G_x(i,j) + G_y(i,j) - G_x(i-1,j);
end
 
for i = 2:imh-1
j = imw;
e = e+1;
A(e,im2var(i,j)) = -4;
A(e,im2var(i-1,j)) = 1;
A(e,im2var(i,j-1)) = 1;
A(e,im2var(i+1,j)) = 1;
 
b(e) = G_x(i,j) - G_x(i-1,j) - G_y(i,j-1);
end
 
e = e+1;
A(e, im2var(1,1)) = -4;
A(e, im2var(2,1)) = 1;
A(e, im2var(1,2)) = 1;
b(e) = G_x(1,1) + G_y(1,1);
 
e = e+1;
A(e, im2var(1,imw)) = -4;
A(e, im2var(1,imw-1)) = 1;
A(e, im2var(2,imw)) = 1;
b(e) = G_x(1,imw) - G_y(1,imw-1);
 
e = e+1;
A(e, im2var(imh,1)) = -4;
A(e, im2var(imh-1,1)) = 1;
A(e, im2var(imh,2)) = 1;
b(e) = G_y(imh,1) - G_x(imh-1,1);
 
e = e+1;
A(e, im2var(imh,imw)) = -4;
A(e, im2var(imh-1,imw)) = 1;
A(e, im2var(imh,imw-1)) = 1;
b(e) = - G_x(imh-1,imw) - G_y(imh,imw-1);
 
 
v = A\b;
 
output = zeros(400,600);
n = 0;
for i = 1:imh
for j = 1:imw
n = n+1;
output(i,j) = v(n);
end
end

7. Convert to the RGB color space

xyz(:,:,2) = exp(imresize(output,10));
bonus_tmp = xyz2rgb(xyz);
imshow(bonus_tmp)
function [output] = weighten(type,im_type)
 
if (im_type == 'jpg')
x = (0:255)/255;
len = 255;
elseif (im_type =='raw')
x = (0:65535)/65535;
len = 65535;
end
if (type == 'unif')
output = ones(len+1,1);
 
elseif (type == 'tent')
for i = 1:len
if x(i) <= 0.5
x(i) = x(i)+0.001;
else
x(i) = 1.001-x(i);
end
end
output = x;
 
elseif (type == 'gauss')
x = exp(-4*(x-0.5).^2/(0.5.^2));
output = x;
end
end
 
function [Z] = image_stack(images,c,s)
length = size(images{1},1);
width = size(images{1},2);
Z = zeros(s,size(images,2));
 
idx_pool = 1:length*width;
sample_set = randsample(idx_pool,s);
 
for i = 1:size(images,2)
tmp = reshape(images{i}(:,:,c),[length*width,1]);
Z(:,i) = tmp(sample_set);
end
end
 
 
function [g, lE] = gsolve(Z,B,l,w)
n = 256;
A = sparse(size(Z,2)*size(Z,1)+n+1,n+size(Z,1));
b = sparse(size(A,1),1);
 
k = 1;
for i = 1:size(Z,1) %픽셀
for j = 1:size(Z,2) %이미지
wij = w(Z(i,j)+1);
A(k,Z(i,j)+1) = wij;
A(k,n+i) = -wij;
b(k,1) = wij * B(j);
k = k+1;
end
end
 
A(k,129) = 1;
k = k+1;
 
for i = 1:n-2
A(k,i) = l*w(i+1);
A(k,i+1) = -2*l*w(i+1);
A(k,i+2) = l*w(i+1);
k = k+1;
end
x = A\b;
g = x(1:n);
lE = x(n+1:size(x,1));
 
end
 
function [G_pyramid] = genpyramid(img,level)
G_pyramid = cell(1,level+1);
G_pyramid{1} = img;
for i = 2:level+1
G_pyramid{i} = G_filter(G_pyramid{i-1});
end
end
 
function [filter_img] = G_filter(img)
g_filter = [1 2 1; 2 4 2; 1 2 1]/16;
temp_img = imfilter(img,g_filter,"replicate","same");
filter_img = zeros(ceil(size(temp_img,1)/2),ceil(size(temp_img,2)/2),size(temp_img,3));
 
ttemp_img = temp_img;
filter_img = ttemp_img(1:2:size(temp_img,1),1:2:size(temp_img,2));
end